home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_9 / rkf45.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  2.0 KB  |  70 lines  |  [MATF/MATL]
  1. function [T,Y] = rkf45(f,a,b,ya,m,tol)
  2. % [T,Y] = rkf45(f,a,b,ya,m)
  3. % Runge-Kutta-Fehlberg solution for y' = f(t,y) with y(a) = ya.
  4. % f  is the function, input.
  5. % a  is the left endpoint, input.
  6. % b  is the right endpoint, input.
  7. % ya is the initial condition, input.
  8. % m  is the initial guess for steps, input.
  9. % T  is the vector of abscissas, output.
  10. % Y  is the vector of ordinates, output.
  11. a2 = 1/4; b2 = 1/4; a3 = 3/8; b3 = 3/32; c3 = 9/32; a4 = 12/13;
  12. b4 = 1932/2197; c4 = -7200/2197; d4 = 7296/2197; a5 = 1;
  13. b5 = 439/216; c5 = -8; d5 = 3680/513; e5 = -845/4104; a6 = 1/2;
  14. b6 = -8/27; c6 = 2; d6 = -3544/2565; e6 = 1859/4104; f6 = -11/40;
  15. r1 = 1/360; r3 = -128/4275; r4 = -2197/75240; r5 = 1/50;
  16. r6 = 2/55; n1 = 25/216; n3 = 1408/2565; n4 = 2197/4104; n5 = -1/5;
  17. big = 1e15;
  18. h = (b-a)/m;
  19. hmin = h/64;
  20. hmax = 64*h;
  21. max1 = 200;
  22. Y(1) = ya;
  23. T(1) = a;
  24. j = 1;
  25. tj = T(1);
  26. br = b - 0.00001*abs(b);
  27. while (T(j)<b),
  28.   if ((T(j)+h)>br), h = b - T(j); end
  29.   tj = T(j);
  30.   yj = Y(j);
  31.   y1 = yj;
  32.   k1 = h*feval(f,tj,y1);
  33.   y2 = yj+b2*k1;                         if big<abs(y2) break, end
  34.   k2 = h*feval(f,tj+a2*h,y2);
  35.   y3 = yj+b3*k1+c3*k2;                   if big<abs(y3) break, end
  36.   k3 = h*feval(f,tj+a3*h,y3);
  37.   y4 = yj+b4*k1+c4*k2+d4*k3;             if big<abs(y4) break, end
  38.   k4 = h*feval(f,tj+a4*h,y4);
  39.   y5 = yj+b5*k1+c5*k2+d5*k3+e5*k4;       if big<abs(y5) break, end
  40.   k5 = h*feval(f,tj+a5*h,y5);
  41.   y6 = yj+b6*k1+c6*k2+d6*k3+e6*k4+f6*k5; if big<abs(y6) break, end
  42.   k6 = h*feval(f,tj+a6*h,y6);
  43.   err = abs(r1*k1+r3*k3+r4*k4+r5*k5+r6*k6);
  44.   ynew = yj+n1*k1+n3*k3+n4*k4+n5*k5;
  45.   if ((err<tol)|(h<2*hmin)),
  46.     Y(j+1) = ynew;
  47.     if ((tj+h)>br),
  48.       T(j+1) = b;
  49.     else
  50.       T(j+1) = tj + h;
  51.     end
  52.     j = j+1;
  53.     tj = T(j);
  54.   end
  55.   if (err==0),
  56.     s = 0;
  57.   else
  58.     s = 0.84*(tol*h/err)^(0.25);
  59.   end
  60.   if ((s<0.75)&(h>2*hmin)), h = h/2; end
  61.   if ((s>1.50)&(2*h<hmax)), h = 2*h; end
  62.   if ((big<abs(Y(j)))|(max1==j)), break, end
  63.   mend = j;
  64.   if (b>T(j)),
  65.     m = j+1;
  66.   else
  67.     m = j;
  68.   end
  69. end
  70.